Spectral dimensions of hierarchical scale-free networks with shortcuts 
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The spectral dimension has been widely used to understand transport properties on regular and 
fractal lattices. Nevertheless, it has been little studied for complex networks such as scale-free and 
small world networks. Here we study the spectral dimension and the return-to-origin probability 
of random walks on hierarchical scale-free networks, which can be either fractals or non-fractals 
depending on the weight of shortcuts. Applying the renormalization group (RG) approach to the 
Gaussian model, we obtain the spectral dimension exactly. While the spectral dimension varies 
between 1 and 2 for the fractal case, it remains at 2, independent of the variation of network 
structure for the non-fractal case. The crossover behavior between the two cases is studied through 
the RG flow analysis. The analytic results are confirmed by simulation results and their implications 
for the architecture of complex systems are discussed. 

PACS numbers: 89. 75. He, 05.40.-a,05.10.Cc 



I. INTRODUCTION 

The problem of random walks (RWs) on complex net- 
works have attracted much attention as a model to study 
diffusion processes on complex systems such as fad or dis- 
ease spreading over social networks, data packets trans- 
port in the Internet, data mining on the web, and so 
on [H-@| . Such importance in theoretical and application 
as pec ts has led to the extensive studies on the problem 
[7H15| , but it still remains unclear how RW properties de- 
pend on network structure. In this paper, we investigate 
the effect of shortcut links on RW motion. Shortcuts may 
be classified into two types: Long-range shortcuts, uti- 
lized in a skeleton or superhighways, and local shortcuts, 
complementary to the skeleton or local roads It 
was recently shown systematically [19( that long-range 
shortcuts can change the type of networks from fractal to 
non-fractal. The fractal (non-fractal) network is the one 
in which the mean separation between two nodes scales 
with system size in a power-law (logarithmic) manner. 
We are particularly interested in how RW properties are 
changed as long-range shortcuts are added and thereby 
the network changes from a fractal to a non-fractal net- 
work. 

The Gaussian model is useful to study RW problems 
analytically [2(|, which we use to understand the RW 
motion on hierarchical scale- free networks [2ll - |24j |. In 
this artificial network, we can control the degree distri- 
bution and the link weight of shortcuts. We consider 
RWs on this hierarchical network, and obtain the exact 
solution of the return-to-origin (RTO) probability by ap- 
plying the renormalization group (RG) approach to the 
Gaussian model. The RTO probability is related to the 
free energy of the Gaussian model and the spectral den- 
sity function of the Laplacian matrix, all of which are 
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characterized by the spectral dimension d s 0, [2(| ■ The 
hierarchical network is recursively organized, so that one 
can obtain analytically the spectral density function by 
using the decimation method of the RG transformation. 
We find that when the network structure is a fractal, the 
spectral dimension d s varies in the range 1 < d s < 2; 
however, it is fixed at d s — 2 for the non-fractal net- 
works. These analytic results are confirmed by numerical 
simulations. Moreover, we perform the numerical simu- 
lations of RWs on various real- world networks and check 
the relationship between the network fractality and the 
spectral dimension. 

The paper is organized as follows. In Section|TTl the hi- 
erarchical network model and the formalism to derive the 
spectral density function of the Laplacian matrix from 
the Gaussian model is briefly introduced along with their 
relation to the RTO probability of RWs. In Section UIT1 
we apply the RG approach to the Gaussian model and 
obtain the spectral density function explicitly. In Sec- 
tion IIV1 numerical results of the RTO probability are 
presented and compared with the analytic solutions. We 
summarize our findings and discuss their implications in 
Section El 



II. MODEL AND FORMALISM 

The network we use in the paper is a modified ver- 
sion of an existing hierarchical network model (2~ir{2~4l j. 
in which shortcuts are assigned weights. We call this 
modified model the weighted flower (WF) network. The 
WF network contains two types of links, type I and II, 
which evolve differently. The WF network is constructed 
iteratively as follows. We start from a single type-I link 
between two nodes, which is the 0-th generation. At the 
next generation, two chains are added between the two 
nodes at both ends of the type-I link. On one chain, u— 1 
nodes are located and thereby u links are placed. On the 
other chain, v — 1 nodes are located and thereby v links 
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FIG. 1. Weighted Flower (WF) network and the renormal- 
ization of the Gaussian model, (a) Construction of the WF 
network with u = 2 and v = 4 from generation n = to 
n = 2. Type-I links (solid line) and type-II links (dotted line) 
have link weight 1 and p, respectively, (b) Renormalization 
of a leaf and a type-II link. A leaf has two chains and one 
type-II link. The upper chain has u — 1 4> variables and the 
lower chain has v — l<f> variables. After a RG transformation, 
each leaf is replaced by a type-I link. A type-II link outside 
the leaves remains after a RG transformation. 



are placed. Those u + v links are assigned as type-I links, 
and the type of the seed link is changed to type II. These 
steps - adding nodes and links and assigning link types - 
are performed at every generation. Every link is assigned 
weight: 1 for each type-I link and p for each type-II link. 
The link weights will be used for the transition rate of 
RWers. The growth of the WF network with u = 2 and 
v = 4 is depicted in Fig.Q](a). We consider only the case 
u > 2 and v > 2 throughout this paper. The type-II links 
play a role as shortcuts in the system and the parameter 
p represents the transition rate across a type-II link. We 
point out that in the previous model of the hierarchical 
network (2l| - [23 |. the parameter p was used as the occu- 
pation probability of the type-II links, but here, it is used 
as the weight of the type-II links. Thus, the two network 
models are equivalent if p = or p = 1. Otherwise, they 
are different. 

The link weights of the WF network at the n-th gen- 
eration are represented in a symmetric matrix form 



W lq 



where 2#> and E<® 



(^)g4 i} , 

(£q) G £i H) 
otherwise, 



(1) 



number of type-I links 
number of type-II links 
number of nodes 
type-I degree 
type-II degree 
node degree 
degree distribution 



W = (« + v y 
(ii) _ (u+vr- 



u-\-v- 



^n — n£-\-l 



k t = (i + p)2 n " 



2p 

"f+2 „ 2p 



p(k) ~ k 7 with 7 = 1 + 







TABLE I. Basic properties of the WF network at generation 
n. Here I is the node index and ne denotes the birth gener- 
ation of node £. The degree of a node is defined here as the 
sum of the weights of connected links. 



links at the n-th generation, respectively, and (£q) repre- 
sents the link connecting node I and q. Basic structural 
properties of the WF networks can be obtained analyt- 
ically [22hl24l | and are summarized in Table HI Here we 
defined the degree kg of node £ as kg = J2 q W q i. The 
distance between two nodes was defined as the length 
of the minimum-cost path with the link cost C£ q given 
by a decreasing function of Wt q satisfying cg q — > oo as 
Wi q —> [25[. While the degree exponent depends on u 
and v but remains the same under the variation of the 
weight p, the mean distance between nodes depends on 
p. When p = 0, the type-II links are not considered in 
determining the minimum-cost paths. Then the mean 
distance D between nodes in the WF networks is related 
to the total number of nodes as \22 241 



D - N 1/d t with d f 



log(w + v) 
log(min{u, v}) 



(2) 



Here df is the fractal dimension (2(| of the WF networks 
with p — 0. When < p < 1, shortcuts may participate 
in the minimum-cost paths and the mean distance is pro- 
portional to the number of nodes logarithmically [22|, [2?J 



D - log N. 



(3) 



Let us consider a RWer located initially at node £q and 
going around in the network with the transition probabil- 
ity Wi q / J2e Wi'q from node q to £. Then the probability 
Pu {t) to find the RWer at node £ after time t is given 
by (J — V)\g Q , where / is the identity matrix and L is 
the Laplacian matrix defined as Li q — Si q — Wi q /k q with 

kq = J2l W tq- 

The RTO probability P (t) = AT -1 E^-P«(*) is deter- 
mined by the eigenvalues fa (i = 0, 1, 2, . . . , N — 1) of the 
Laplacian matrix L as 

N-l 



P (t) = iTr (/-£)' = I £(1- Mi ) 



denote the set of type-I and type-II 



d/i />(//) e-^ + (-l)* / djipp-Me-t*, 
o Jo 

(4) 

where we introduced the spectral density function p(fj,) = 
A^ 1 ^^q 1 S(fi — fa) in the thermodynamic limit N — >• 
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oo [28{. For large t, the term (1 — pif in Eq. (|4]) is 
dominant as pi — > and pi — > 2. The contribution 
from fa ~ 2, however, cannot be larger than that from 
|Uj ~ since P (t) > 0, but can cause an oscillatory be- 
havior in P {t) depending on whether t is even or odd. 
For example, in d-dimensional regular lattices, the eigen- 
values are symmetrically distributed about p = 1, and 
p(p) = p(2 - p) ~ p d l 2 ~ l for small p. Then, it fol- 
lows from Eq. §4§ that P {t) ~ t~ d / 2 when t is even and 
P Q (<) = when t is odd. 

For fractal lattices, the spectral density function and 
the RTO probability are characterized by the spectral 
dimension d s as @: 



p(p) ~ p?'' 2 - 1 for p—> Q, and 
Po(t) ~ i~ ds/2 for t -> oo. 



(5) 



We use these formulae in further discussion later. 

The mean first passage time T between two nodes is 
related to the eigenvalues of the Laplacian matrix via 
T = P^ 1 Q, which scales with system size N as T ~ 
N 2/d B for rf s < 2 and T - iV for d s > 2 [ll|. For the 
fractal WF network with p = in the n-th generation, 
it was shown that T ~ (uv) n ~ 7V lo s(™)/ i°s(«+*) baged 
on the scaling behavior of the mean first passage time 
between two hub nodes 23]. Thus, it follows that d s = 
2 log(u + v)/ log(uu) for the case p = 0. On the other 
hand, the mean first passage time for a non-fractal WF 
network with u = 1, v = 2 and p = was calculated, 
which scales as T ~ 3™ ~ N [Ul^. Thus, the spectral 
dimension is d s = 2 in this case. These results suggest 
that the spectral dimension can be different depending on 
whether a network is fractal or non-fractal. Motivated by 
these previous works, we study the spectral dimension of 
the hierarchical network as it crosses over from a fractal 
to a non-fractal by varying the parameter p in the range 
< p < 1. 

To derive the spectral density function of the matrix L, 
we use the Gaussian model of which the partition func- 
tion is given as 



Z(p) 



T>4> exp 



: ,22(f>eHt q (p)<j> q 



N 



(^) 

det.ff' 



r N 



(6) 



and 



where N is the number of nodes, T>4> — ]^L=i ' 
H = pi — L. L is a similarity-transformed Laplacian 
matrix given in Appendix [A] p is a parameter intro- 
duced in the Gaussian model. The spectral density 
function p(p) is derived from the free energy f(p) = 
- limAr^oo N^ 1 log Z(p) as 



( ^ 2 T d f 

TT Op 



(7) 



where p is assumed to possess a positive infinitesimal 
imaginary part. 



III. RENORMALIZATION GROUP APPROACH 
TO GAUSSIAN MODEL 

A. Fixed points and RG flow 

The parameter p in Eq. ([5]) is renormalized differently 
for the two types of links in the WF networks. Thus, we 
first rescale <pi as 4>i / sfkt in Eq. © and introduce distinct 
parameters p\ and P2 for the respective link types. The 
Hamiltonian then takes the form H? q = [(pi — l)k^ + 

(P2~ ]Siq + Wiq, where Wt q is given in Eq. ([TJ. The 
partition function in the n-th generation is written as 

Z n {pi,P2,p) 

=jv4> n iPifah) n 4 n) (^,^). 

where 

i$(M q ) = exp (-t(l - pi)($ + <f> 2 q ) + 2*0/0,) , 
and 

4 H) (^> <M = ex P (~*(! ~ P>2)P{4>1 + 4> 2 q ) + 2ip<j>t<l>q) ■ 

(8) 

Note that p\= pi = p initially. The RG transformations 
for the shortcut weight p, and the parameters pi, and p% 
are carried out to obtain the spectral dimension d s . As we 
show later, p 2 is invariant under the RG transformation 
and we investigate the fixed points and the RG flow in 
the (pi , p) plane. If there is only one fixed point at p = 0, 
then the spectral dimension d s will be the same regardless 
of p. On the other hand, more than one fixed points will 
indicate a crossover or a transition of d s depending on p. 
The WF network in the n-th generation contains L^ > _ 1 

leaves and L^ I } 1 links, and the partition function may be 
written as 

z n =fv4>' n ZtetWtA) n &(M q )- 

(1\ /TT^ 



(i'q')eE^ 



Here, each leaf contains u+v nodes, among which u+v—2 
are the youngest nodes, and the two remaining old nodes 
are connected by a type-II link as shown in Fig. [lib). 
Integrating out the 4> variables on the youngest nodes, we 

find that Z\ ea s((j)' 1 , 4>' 2 ) takes a similar form to ^liO^i , <p' 2 ) 
as 

Z leaf (^,^) = («/2)("+^ 2 )/ 2 G(u 1 )- 1/2 x 
exp 



-% {2(1 - Ml ) + P {i - p 2 ) - /n^)} (4? + 4>' 2 2 ) 

+2i{p + fta(Mi)}M, (9) 
where 

c(£, a, b) = 1 — a — cos(^7r/6), 
G(pi) 



u—1 v—1 

Y[ c(i,p!,u) Y[ c{q,pi,v), 

9=1 



4 



u — 1 .9//? / \ i u ~ 1 -2/ /\ 

, . 1 r-s Sin (lir/u) 1 ^-^ sm (qTr/v) 

U fcl C (*>Ml> U J V ' ~y C (<7>Ml,«) 



and 



l^. 1 (-I)*-* sin 2 (^) 



D-l 



lyi (-l^sin^gTr/u) 



■7=1 



c(q,m,v) 



(10) 



The derivation of Eq. (|9]) and the properties of 
G(fj.i),hi(ni), and ft.2(Mi) f° r small /ii are given in Ap- 
pendices iBland ICl In order for Zi oa f 02 ) m Eq. (|9]) 
to take the same functional form as 0' 2 ), all the 

4> variables should be rescaled as cf>' / \fp + h^pi) so that 
the coefhcient of 4>'i4>2 be unity. With the introduction 
of n[ , \i! 2 , and p' given by 

, 2(l- / u 1 )+p(l-/z 2 )-/ii(Mi) 
1 - Mi = 



p + h 2 (ni) 



M 2 = M2, 



(11) 



J3+/l 2 (/Jl)' 

the partition function Z ra is related to Z n ^\ as 

Z n (pufi2,p) = exp[-Af n g(/ji,p)]Z n _i[/j' 1 ,/j2 / ,p'], (12) 
where the function g(n\,p) is defined as 

ff(Mi.P) = log[G(Mi)] + -^-log[p + /i 2 (Mi)] ■ 

" ~" (13) 

As this RG transformation is repeated, the parameters /xi 
and p are renormalized successively according to Eq. ([lip . 
The parameter /j 2 remains at its initial value fi. On the 
other hand, the parameter \x\ is renormalized in a non- 
trivial way. 

For small /ii, one finds two fixed points in the (ni,p) 
plane located at 

(jj,l,p*) = (0,0) and (0, 1 - u' 1 - iT 1 ). (14) 

These two fixed points meet when (u, v) — (2, 2). Let us 
first consider the case (u, v) ^ (2, 2). Near the fixed point 
(fj,* = 0,p* =0), the RG equation, (fTTl) . is linearized as 



/i'j =uv(j,i and p 



(/(> 



U + V 



p. 



(15) 



and therefore, after r RG transformations, /ii and p take 
the following values: 



Hi = (uv) T fi and p^ — 



u + v 



(16) 



Since uv and uv/(u + v) are not smaller than 1 for u > 2 
and v > 2, the fixed point (0,0) is unstable. Around 



the other fixed point (fi* = 0,p* = 1 
linearized RG equation is written as 



the 




R 



Mi 

P~P* 



M 



1 — u — v 



R 



u + v 

(uv — u — v)(uv — l)(u+v) u+v 
3(uv) 2 uv 



(17) 



The matrix R has two eigenvalues u + v and (u + v)/(uv) 
and the parameters after r RG transformations are given 
by 

(u + v)(uv — 1) 



(t) 

Mi 



M 



uv(u + v — 1) 
(u + u)(fiv — l)(uu — u — t>) 



+ (p-p*) 



3(uu) 2 (u + u 
u + v 



{(u + v) T - 1} + 1 

{i-(«+t,n 

(18) 



1) 



If (u, v) — (2, 2), there is only one fixed point at (fi* = 
0,p* = 0) for small The renormalized value p^ can 
be calculated from Eq. (fTTj) as 



(19) 



Note that p( T ) is not as small as /j. Inserting Eq. flT9")) into 
Eq. dlip . a recursive relation for fi\ is obtained, which can 
be solved to give 



(t) 

Mi 



(1 + 3p- 1 )4 r - 1 



: M- (20) 

3(r + p i ) 

To understand the renormalization of /ii and p, which 
depends on their initial values as well as the the net- 
work model parameters u and v, we consider the follow- 
ing three cases: 

(I) P = 0, 

(II) p^0 and (u,v)?(2,2), 

(III) p^0 and (u,v) = (2,2). 

(I) If p = initially, p remains at zero and /ii in- 
creases according to Eq. (|T6"|) . (II) If p > initially 
and (u, w) ^ (2,2), both /ii and p increase by Eq. (fTo| 
and the RG flow runs outward from (0, 0) along the 
curve p ~ /ij log («+■")/ log (™) ag obtained from Eq. (|16p . 
The RG flow then runs towards the other fixed point 
(0, 1 — it -1 — u _1 ), driven by one eigenvalue (u+v) / (uv) < 
1 of R in Eq. (TT7|) . Once p gets sufficiently close to p*, 
the RG flow is bent outward along the line p — p* — 
—fil(uv — (u + v))/(3uv), parallel to the other eigenvec- 
tor of R (see Fig. [1(a)). (Ill) When (u,v) = (2,2), the 
two fixed points coincide each other. The RG flow along 
the p axis runs quite slowly as the corresponding eigen- 
value of R is 1. If p = 0, [i\ and p are renormalized in 
the same way as in (I). If p ^ 0, the RG flow runs along 
the curve fi\ ~ p4 1 / p as obtained from Eq. (|2U)) . which is 
depicted in Fig. &b). 
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(a) 
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v-u-v) JJ.j 



(0,1- 
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(b) 
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FIG. 2. RG flow in the (pi,p) parameter plane, (a) RG flow 
for [u, v) 7^ (2,2). There are two fixed points in the small-^ii 
region: (0,0) and (0, 1 — u _1 — If p = initially, the 

flow is outward along the /ii axis. When p / initially, it is 
first attracted towards (0, 1 — u^ 1 — v^ 1 ) and then repelled 
along the line p = 1 — (it + v)/(uv) — (j,i(uv — u — v)/(3uv). 
(b) RG flow for (u,v) = (2,2). There is a fixed point at the 
origin and the flow follows a curve p,\ ~ pA 1 ^. 



B. Spectral density function and spectral 
dimension 

The free energy per node in the thermodynamic limit 
f(n) = -lim N n 1 log Z n (ii /i, p) can be obtained by 
using g(pi,p) in Eq. (fT3)) as [20| 



m = E 



i (« + «) T 



(21) 



where we used N n -\/N n ~ {u + v) 1 for large n. Using 
the expansion of the functions G(/ii) and /la(Mi) given in 
Appendix [Bl we can obtain the expansion of g^i^.p^) 
up to the first order in and p^ — p* as 



2(u + v)(u + v - 2) 



log(G(0)) + log 1 




1 



2(u + v) 



[\og[G(0){p* +U- 1 +V- 1 }] 



log 1 + 



o(r) 



P" 



,,(t) (uv-1){u+v) (t)' 
Ml 3uv Ml 



Mel 



,(22) 



where we used L n _i/N n ~ (u + u— l)/((n+u)('U+u — 2)) 
and introduced a constant 



Pel 



(23) 



The renormalized values of /ii and p obtained above for 
each case of (I), (II), and (III) can be used in Eqs. (|2"Tj) 
and ([22]) to evaluate the free energy /(/x). 

For the case (I), the shortcuts have no weight and the 
network is a fractal as long as u > 2 and v > 2. The pa- 
rameter p remains at zero under the RG transformation 
and fix is renormalized by Eq. (|16[) . Inserting Eq. (|16[) in 
Eqs. dlT|) and (JUJ) with p( T ) = p* = 0, we obtain the free 
energy and then the spectral density function via Eq. (JT]) , 
which is calculated as 



M(i)(m) 



7T 



S(z-n al /n) 
(u+v)(u+v-2) 



5(z-fj, a2 /n) 
u+v 



PL \0g(uv)z l °s( u + v )/ ^S(t- 



(24) 



where 



and 



Mc2 = 



v 2 — — 1 ' 



K m 



((«• 



log(u + ^) 
-1 l og (m>) 
Mcl 



M c2 



(u + v) \og(uv) 



Here we used the relation logx = log \x\+iir6(— x) for real 
x, where 0(x) is the Heaviside-step function. Comparing 
Eq. ([5]) with ([24]) . we find that the spectral dimension d s 
for the case (I) is 



d {i) = 2 logpu + v) 
log(ui;) 



(25) 



This is consistent with the previous result obtained for 
the mean first passage time studied in Ref. [23]. The 
spectral dimension in Eq. (I2~5j) is reduced to 2 when 
U = v = 2, and ranges between 1 and 2 when u and 
v differ from 2. Note that it is different from the fractal 
dimension given in Eq. 

For the case (II), the RG flow starting from a point 
with non-zero values of pL\ and p first approaches the 
fixed point (0,1 — u^ 1 — v^ 1 ), and then it is repelled 
in the direction of one eigenvector of the matrix R in 
Eq. (fTTl) as shown in Fig.[^a). Using Eq. (IT51) and setting 
p* = 1 - ir 1 - v' 1 in Eqs. ([21]) and ([22]), we obtain the 
spectral density function 



M(ii)(m) 



2 

T = ^ 



9 4t^ t) ,pW) 



7T 



(m + w) 7 



6 



(u+v)(u+v-2) 



si*-*?) 

(u+v) 



' 1 



Zjjj log(u + v) 



(26) 



where 



Pc3 = 



(it + v)(uv — 1) 
(uv)(u + v — 1) 



3mw(u + u — 1) 



and 



(iti> — l)(it + t>)(u 2 + v 2 — u — v — 1) ' 
(uv — l)(u 2 + v 2 — u — v) 



MID 



3uv (u + v — 2) log(u + v) 



Thus, the spectral density function for the case (II) is a 
constant for small /i, indicating that 



di n > - 2. 



(27) 



This means that the spectral dimension is changed by 
the addition of shortcuts. Remarkably d s remains at 2 
regardless of u or v. 

For the case (III), there is only one fixed point and 
p( T ) remains to be 0(1). Inserting Eqs. (fT9|) and ([20 



in Eqs. (|2"Tj) and (|2"2l with p* — 0, the spectral density 
function for u = v — 2 is calculated as 



2 

P(iii)(m) - -Im, 



oo 



(«• 



log 4 



1 



Mci 3(1 +plog 4 z) 



1 



4 

P + 3 



1 3(1 
1 



f>log 4 z) 



3 



61og4 1+plog 4 (l/ A i)' 



(28) 



We remark that Eq. (|28[) with p = is reduced to 
Eq. On the other hand, for p > 0, the spec- 

tral density function takes a logarithmic form, p(fJ>) ~ 
(log(l//i)) _1 for /it — > 0. In the scaling regime where 
/iexp(p _1 ) is finite, the spectral density function displays 
the following crossover behavior: 



P(iii)(m) 



nr^i — ttt^ for /ie p 1 < 1, 
6 P+ io g4 (i/M) ^ » ^ ^ > (29) 



6 log 4 



for /j,e p > 1. 



IV. THE RETURN-TO-ORIGIN PROBABILITY 

In this section, we derive the RTO probability P (t) 
of RWs from the spectral density function p(fi) obtained 
in the previous section. The behavior of p(/j.) around 
fi = 2, which should be known for Eq. ((3]), can be un- 
derstood by the symmetry property of the RG equation 



(fTTT) . First let us consider two cases, (i) When both u and 
v are odd, G(iii) and /^(mi) are symmetric and h\{pi) 
are anti-symmetric with respect to fix = 1. Thus, the 
variable p,i = 2 — /xi and p are renormalized in the same 
manner as Eqs. (|T6)) and (p~8|) around the fixed points 
{fi\ = 2,p* = 0) and {fi\ = 2,p* = 1 - u~ x - v~ v ), 
respectively. This leads to p(2 — /z) = /»(//) ■ (ii) When 
p = and both it and v are even, G{pi) are symmetric 
and hi(pi) and /i2(pi) are anti-symmetric with respect 
to /iti. Interestingly, the parameter fix initially close to 
2 is renormalized by one RG transformation to a value 
near zero by Eq. (fTTj) , and then follows the same RG flow 
as given in Eqs. (fl"6|) and (fl"8]) . Therefore, the spectral 
density function is symmetric about (J, = 1. In these two 
cases (i) and (ii), it follows that 



P (t)~(l + (-l) ( ) / dwMe 



,-fj.t 



(30) 



For other cases than (i) or (ii), there is no fixed point on 
the line p,\ — 2 and thus the contribution of p{p) around 
p, = 2 to P (t) for large t is subdominant, leading to 



Po(t) 



dpp(p)e 



(31) 



Taking account of different behaviors of p(p) for small 
p, in the cases (I), (II), and (III), we obtain from Eq. (|U) 

f r(d«/2) ^ dlI)/2 (i), 

P (i)~ P J ^(n)*- 1 (II), (32) 

I ^^(i+pfog^r 1 (ni), 

where r(x) is the Gamma function and the spectral di- 
mension dffl in the case (I) is given in Eq. (|2"5j) . Note 
that for the case (III), a logarithmic term appears. As 
discussed above, the coefficient p a is 2 when both u and 
v are odd or when p = and both u and v are even. 
Otherwise, it is 1. 

We performed numerical simulations of RWs on the 
WF networks. In the simulation, 10 6 independent RWers 
moved around on the networks with the transition prob- 
ability Wi q /k q in Eq. (fTJ. The RTO probability P a (t) is 
measured as the fraction of the walkers who are found 
in their own initial positions at t. The measured RTO 
probabilities for various values of p and (u,v) are shown 
in Fig. O 

The spectral dimension is obtained from the power-law 
decay of P (t). When p = 0, the measured values of d s 
for various u and v are in good agreement with Eq. (|25p 
(see Fig. [3] (a)). On the other hand, when p = 1, the 
numerical values of d s are somewhat larger than the the- 
oretical value 2, weakly dependent on u and v ranging 
2 ~ 2.2 as shown in Fig. [3] (b). Given the logarithmic 
term for u = v = 2 in Eq. (|3"2")l . this slight deviation as 
(u,v) approaches (2,2) can be a finite-size effect. The 
faster decay of P (t) for p = 1 than p = is caused by 
shortcuts. Other structural factors such as the degree 
distribution, however, do not affect the spectral dimen- 
sion, unlike the case p = 0. 
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FIG. 3. Return-to-origin probability P (t) on WF networks, (a) Plots of P (t) on the fractal WF networks (p = 0). Simulation 
results for (u, v) = (2,3) and (u,v) = (5,5) are shown and the lines represent the theoretical prediction in Eq. (|32[1 . (inset) 
Plot of the spectral dimension d a as a function of 21og(u + v)/log(uv). The value of d a were estimated from P (t) for 
(u,v) = (2, 3), (2, 4), (2, 6), (3, 3), (3, 4), (3, 5), (3, 6), (4, 4), (4, 5), (4, 6), and (5,5) using the relation P (t) ~ t~ d ° /2 . The line 
represents the analytic prediction d s — 21og(u + «)/log(wv). (b) Plots of P a (t) on the non-fractal WF networks (p — 1) for the 
same values of u and v. The lines represent the theoretical predictions in Eq. I|32|l. (inset) Plot of the spectral dimension d a as 
a function of 2 log(« + v)/ \og(uv). The line represents the analytic prediction d s = 2. (c) Plots of P {t) for (u,v) = (3,5) and 
various values of p. The slopes of the two lines are — log 8/ log 15 and — 1, corresponding to the analytic prediction for d s /2 for 
p = and p = 1 for (u, v) — (3, 5). (inset) Plot of the spectral dimension d s as a function of p. The results obtained for the 
5th and 6th generation are presented, (d) Plots of (1 +plog 4 t)P (t) / P ,th(X) versus t for (u,v) = (2, 2) and p = 0, 0.2, 0.5, 0.8, 
and 1, where P ,th(l) = Po(p + 3)/(6 log 4). All the data points collapse to t' 1 as predicted by Eq. (|32[) . 



The case < p < 1 is also interesting. For (u, v) ^ 
(2,2), we obtain that d s increases rapidly with p as 
shown in the case (u,v) — (3,5) in Fig[3] (c). The 
curve looks more concave as TV increases. Thus we 
expect that d s is independent of p in the thermody- 
namic limit as long as p > 0. For (u,v) — (2,2), the 
RTO probability is expected to have a logarithmic term 
as P (t) ~ P . th (l)^ 1 (l +p\og 4 t)- 1 with P , th (l) = 
Po {p + 3)/(61og4). We plot (1 + plo gi t)P {t)/P , th {l) 
versus t for various values of p in Fig.[3Jd), in which the 
data are well collapsed onto i -1 as predicted in Eq. (|32| . 



V. SUMMARY AND DISCUSSION 

In this paper, we have obtained the exact solution for 
the RTO probability of RWs and the spectral density 
function of the Laplacian matrix in hierarchical scale- free 
networks by applying the RG approach to the Gaussian 



model. Particularly, we monitored the spectral dimen- 
sion as we controlled the weight of the shortcuts. When 
the shortcut weight is zero, the network is a fractal, and 
its spectral dimension varies in the range (1,2]. It de- 
pends on the parameters controlling network structure. 
On the other hand, when the weight is larger than zero, 
the spectral dimension is 2, which is robust under struc- 
tural variations. This result demonstrates that shortcuts 
strongly affect the RW motion even when the shortcut 
weight is small. 

It was shown that the organization of shortcuts in 
a network was essential in the classification of network 
structures [Hj]. Let p(r) denote the probability of find- 
ing a shortcut between a pair of nodes that would be 
apart by distance r. If the decay of p(r) is slower than 
r~ 2d f with df being the fractal dimension of the network, 
the network is renormalized to the complete graph. Oth- 
erwise, the network is renormalized to the corresponding 
fractal network. In the non-fractal WF networks, we ob- 
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tain that p(r) ~ r~ a , where a/df is in the range between 
1.6 and 2. Thus, the hierarchical networks with shortcuts 
belong to the former case, which can be related to our 
finding that the spectral dimension d s is robust against 
the variation of u or v in the p > regime. 

We performed the numerical simulations of RWs on 
the following real- world complex networks: the protein 
interactions network of the yeast [30| . the human pro- 
tein interaction network 1311 . the coauthorship network 
of the cond-mat archive 13211 . and the Internet at the 
autonomous system level [33| . Among them, the two 
protein interaction networks are known as fractals and 
both the coauthorship network and the Internet are non- 
fractals (34)]. The spectral dimensions obtained from the 
power-law decay of the RTO probabilities turned out to 
be larger for non-fractal networks than fractal networks. 
The numerical values of d s were approximately 1.4 and 
2 for the yeast and the human protein interaction net- 
work, respectively, and 4 and 5.1 for the coauthorship 
network and the Internet, respectively. The probability 
R that a RWer will ever return to the origin is related to 
the RTO probability as R = 1 - 1/E^o and thus 
RWs are recurrent (transient) if d s < 2 (d s > 2) [H-Q. In 
this sense, the structure of the hierarchical networks with 
shortcuts studied in this work, having d s = 2, is marginal 
in the RW motion; RWers can never return to the origin 
if the network structure is more entangled. It would be 
interesting to see the case of the Internet. Since the spec- 
tral dimension is larger than 2 for the Internet, packets 
roaming on the Internet may not return to its origin or a 
destination if they follow the RW motion. Therefore the 
studied hierarchical network could be a good sample for 
the construction of the Internet. 
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Appendix A: Other Laplacians 

The Laplacian matrix L, defined as L? q — S£ q — Wi q /k q , 
can be brought into a symmetric form by a similarity 



transformation L = SLS with St, 



,-1/2; 



■-Wla 



(Al) 



The eigenvalues of L and L are identical. Another Lapla- 
cian matrix L' lq = kgSi q — Wg q [35[ has different eigen- 
values from L or L. We consider L in Eq. (|Al|) for our 
analysis of the R W p roblem but the main results hold 
also for L and L' 1361. 



Appendix B: Derivation of Eq. ([9]) 

In this section, we derive the partition function 
Zieaf^i,^) in Eq. © for a leaf in Fig. ffl (b). Using 
Eq. ([8]l. we can represent Zi ea { ((/)[, <p' 2 ) as 

/u— 1 V — 1 
n^ up) n< o) 
n 1 r . 1 



q=l 



u—1 v-2 

A 1 )/J u p~> J u p~)\T1 JD/aHo) Jio) 



n 



■n Wt 



9=1 



x4Vi^i l0) )^-W' 2 ) 



u— 1 v— 1 

(up) 



4 n) (0'i,^) / nM up) iP4 io) 

q=l 



exp 



u—l v—1 

\q=l £,q=i 



E 



+2^(^ p) + <^ o) ) + 20,(^1 + <e\) 



(Bl) 



Here we introduced an (u— 1) x (u— 1) Hamiltonian J$"( u p) 
and a (v — 1) x (u — l) Hamiltonian H^°\ both of which 
take a form 



H 



( H Q -1 \ 

-1 Ho -1 

-I Ho -I ••• 



V 



(B2) 



-1 Ho J 



with Ho = 2(1 — Hi). The (n — 1) x (n — 1) matrices in 
this form have n — 1 eigenvalues Xe and the corresponding 
eigenvectors with I = 1, 2, . . . , n — 1 [37j 



A/ = i?n — 2 cos 



(tt 



— I sin 



, sin 



&r\ /2£tt 
— I , sin ( 

(n - l)£n 



(B3) 



Therefore H^ and #( l0 ) are diagonalized with bases 



I(up) 



1(1°) _ 



u-1 



7 sin 
u z — ' \ u 

9=1 x 

-u-l 



and 



O— 1 v 7 



respectively, and Z\ ea f (<^ 3 </> 2 ) is rewritten as 

Zleaf = exp [-<(1 - /i2)p(0? + ^) 
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/u—X v — 1 

n^ up) n^ io)ex p[ 

1=1 q=l 
1=1 q=l 



u * — ' \ u 

i=i v 



q=l 



-2i(l- M iM 2 +^ 2 )] 



(B5) 



which is calculated to give Eq. ©. Note that c(£, a, b) in 
Eq. (|10p is half the eigenvalue A# in Eq. (|B3[) with n = b 
and i?o = 2(1 - a). 



Appendix C: Expansion of G((j,i),hx(fJ.i), and faOui) 
in Eq. ([TDD 

The small- /m behaviors of G(jJ.i),hi(jii), and /^(mi) 
are used to linearize Eq. (jlip and to find the small-/ii 
behaviors of g(pitP) in Eq. (|13|) . Expanding these func- 
tions for small fj,% up to the first order, we obtain 

1-f {Si(u) + Si(v)}] +0( M 2 ), (CI) 



= G(0) 

u - 1 . u - 1 

AM. 0«l) = 



+ Mi x 

it t> 

Si(«)-(t*-l) , 5i(«)-(«-l) 



/l2(Ml) 



1 1 



! + (-!)" 



&(«) 



(C2) 



where G(0) = n"=i C 1 ~ cos(&r/«)) n^iC 1 ~ cos{qn/v)) 
and we introduced two sums Si(n) and S^iri) defined as 



n— i | 

Si(n) = > — „ — ; — • 



(C3) 



To compute Si(n) and ^(n), we use the following iden- 
tity 



_L = J_ ~ 1 

sin 2 (7ra) tt 2 . ^ (fc - 



k=-c 



(fc - a) 2 



(C4) 



which can be derived by considering the following integral 
on a circle the radius R of which goes to infinity: 



, cot(7rz) 

dZ—y T7T = I 



E 



i 



(z-a) 2 } ^ (n-a) 2 sin 2 (7ra) 

^ n= — oo v ' 

Using Eq. (|C4|) , the sum Si(n) is evaluated as 



= 0. 
(C5) 



n-l 



5*i (n) = > — ^ — : — 
fe-n 2 (^) 



71 — 1 OG 

X, 



i 



i 



j=lfc=-oo \(2n ( X 2n fc ) 

-fl^Vfi 1 A _ 2(n 2 - 1) 



7T ' — ' \ S 

s— 1 x 



(C6) 



Similarly, the sum S2(n) is given by 



H -n 2 (f) 

^ t2 (A 



(-1) 



(j-2nk) 2 {-j -2n{k-l)f 



j = l k=0 

An 2 v^oo ( — 1)° 

TT 2 2-~is = l "s 2 r (sn) 
4n^°° (-1) 3+1 



4ra v^oo 

n 2 - 1 l + (-l) 7 



7 = 3 n is even, 
2 — i 

- = - n is odd, 



3 



(C7) 



Inserting Eqs. (JCS) and (JC7J into Eq. (|C2]) . we find 
that 



u 2 + v 2 - 2 



G( Ml )~G(0) f 1 — a*i 

~ 2 hfii x 

uv 

/2u 2 - 3u+ 1 2v 2 -3v + l 

( 



3w 



3r 



M + U / U — 1 IT — 1 

«2 (a*i) — 1" Mi 



3« 



3i> 



(C8) 
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